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R  ADC  MULTI-DIME  NSIONAL  SIGNAL-PROCESSING 
RESEARCH  PROGRAM 


1.  INTRODUCTION  AND  SUMMARY 

The  Lincoln  Laboratory  Multi-Dimensional  Signal- Processing  Research  Program  was 
initiated  in  FY  80  as  a  research  effort  directed  toward  the  development  and  understanding  of  the 
theory  of  digital  processing  of  multi-dimensional  signals  and  its  application  to  real-time  image 
processing  and  analysis.  Specific  areas  of  interest  include  advanced  filter  and  processor  archi¬ 
tectures  for  image  enhancement,  and  the  modeling  of  image  data  for  segmentation,  classifica¬ 
tion,  and  detection.  This  report  discusses  technical  accomplishments  over  the  last  six  months 
in  these  two  areas. 

Segmentation  and  classification  are  important  aspects  of  the  automated  analysis  of  aerial 
reconnaissance  imagery.  Our  work  on  segmentation  and  classification,  discussed  in  Sec.  2  of 
this  report,  has  focused  on  the  development  of  statistical  models  and  related  algorithms  for 
segmentation/classification.  The  models,  which  are  based  on  white-noise-driven  linear  filters, 
permit  development  of  the  joint  probability  density  function  oi  '  kelihood  function  for  the  image. 

With  an  expression  for  the  likelihood  function  in  hand,  segmentation  can  be  regarded  as  an 
estimation  problem  for  the  image  regions.  Maximum  likelihood  estimation  leads  to  a  simple 
segmentation  algorithm  but  one  that  yields  generally  unsatisfactory  results.  The  modeling  of 
region  transition  statistics  as  a  Markov  chain  allows  formulation  of  the  maximum  ji  posteriori 
(MAP)  estimation  problem;  a  solution  of  the  MAP  equations  can  be  obtained  through  an  iterative 
procedure.  The  MAP  approach  was  found  to  produce  much  more  satisfactory  results.  Examples 
of  image  segmentation  are  shown  for  texture  images  and  for  aerial  photographs  of  rural  areas. 

The  discussion  in  Sec.  3  on  the  enhancement  and  filtering  of  aerial  photographs  represents 
a  coalescence  of  work  in  two  different  areas,  the  development  of  an  iterative  implementation  for 
2- D  non-causal  rational  filters  and  the  investigation  of  techniques  for  restoring  images  that  have 
been  blurred  by  a  known  blurring  function.  We  have  found  that  an  iterative  restoration  algorithm, 
which  is  related  to  the  iterative  implementation  developed  previously,1  works  well  in  removing 
the  effects  of  blurring. 

2.  IMAGE  MODELING  FOR  SEGMENTATION  AND  CLASSIFICATION 

2.1  Introduction 

The  problem  of  detecting  objects  in  structured  backgrounds,  such  as  terrain,  is  of  interest 
in  many  applications.  Since  the  image  in  which  the  detection  is  to  take  place  generally  consists 
of  many  different  types  of  background  regions  with  irregular  boundaries,  segmentation  of  the 
background  into  regions  with  homogeneous  structure  can  be  an  important  prerequisite  to  detec¬ 
tion.  Further,  the  classification  of  background  regions  is  of  interest  since  it  may  reduce  the 
number  of  areas  that  need  to  be  considered  for  detection  of  a  particular  type  of  object. 

During  this  reporting  period,  work  has  focused  on  the  development  of  statistical  models  for 
natural  terrain  and  on  the  generation  of  segmentation  and  classification  algorithms  based  on 
these  models.  The  models  are  derived  using  stochastic  filtering  concepts  and  are  described  in 
the  following  subsection;  the  segmentation  and  classification  algorithms  are  described  in  the 
next  two  subsections.  Examples  of  segmentation  results  for  texture  and  terrain  images  are 
shown  and  discussed. 
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2.2  Modeling  of  Texture  and  Terrain  Images 

Images  of  interest  here  are  assumed  to  consist  of  a  number  of  connected  regions  each  of 
which  contains  a  texture  or  terrain  of  a  fixed  type.  A  model  for  the  image  can  be  developed  as 
follows:  Within  a  given  region,  the  image  is  represented  by  the  output  of  spatial  linear  filter 
driven  by  white  noise  (see  Fig.  1).  If  the  probability  density  function  for  the  white  noise  is  known. 


Fig.  1.  Model  for  image  consisting  of  two  region  types. 

then  the  multivariate  probability  density  function  for  the  set  of  all  points  Fj  within  the  region  R. 
can  be  computed.  Denoting  this  quantity  by  wRere  11  *s  tRe  texture  type,  we  can  write 

the  multivariate  probability  density  for  the  image,  conditioned  on  a  set  of  regions  Rj, . . ,  ,  R^  as 

N 

P(F  |  Rj , . . .  ,  Rn)  =  n  Pk.(FiiRl>  • 

i=l  1 

If  a  causal  (quarter-plane)  recursive  filter  and  Gaussian  white  noise  are  used  in  the  model, 
then  the  expression  for  minus  twice  the  log  of  the  density  function  is  approximately  ^ 

(E  (n.m))2 
V  Ki 

->  2 
n,rm  Hj  CTkj 

where  n  and  m  are  coordinates  of  points  in  the  region,  a,2  is  the  white  noise  variance  and 

Ki 

E^ln.m)  is  given  by 

I  J 

Ek(n,  m)^^  ^  (F(n  —  i,m  —  j)  -  MfcJ  .  (3) 

i-0  j=0 

(k) 

Here  M,  is  the  mean  of  F,  a  -  1,  and  the  remaining  a.  are  the  coefficients  of  the  recursive 
k  oo  >J 

filter.  Under  these  conditions,  Ek  can  be  interpreted  as  the  error  in  linear  prediction  of  the 
image  using  the  filter  of  class  k.  Equation  (1)  and  its  explicit  realization  |Eq.  (2))  are  the  key 
to  the  segmentation  and  classification  algorithms  described  below.  Generalizations  and  further 
details  of  the  modeling  are  given  in  Ref.  2. 

2.3  Segmentation  as  an  Estimation  Problem 

Given  an  image  F,  we  ran  regard  Eq.  (i)  as  a  likelihood  function  in  the  parameters  N  and 
Rj,  A  maximization  of  Eq.  (1)  with  respect  to  these  parameters  yields  a  maximum  likelihood  (Ml.) 

t  The  expression  neglects  boundary  effects  (see  Ref.  2). 
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estimate  for  the  regions.  In  view  of  Fq.  (2),  maximization  of  Eq.  (1)  is  equivalent  to  assigning 
each  point  (n,  m)  in  the  image  to  a  region  (of  type  m)  such  that  the  term  in  brackets  is  mini¬ 
mized.  Specifically,  for  the  case  of  two  region  types,  we  are  led  to  the  condition 

K^jn,  m)  >  0  E^n.m) 

- 2 -  Mn%  - - 2 -  +  lnai  •  (4) 

a  1  u,*" 


The  double  inequality  sign  indic  ates  that  a  point  (n,m)  is  assigned  to  region  U  if  the  left-hand 
expression  is  less  than  the  right-hand  one  and  assigned  to  region  1  if  the  opposite  is  true. 

Since  the  Ml.  estimate  assigns  points  to  texture  types  without  regard  to  the  assignments  of 
adjacent  points,  the  procedure  tends  to  produce  a  poor  estimate  for  the  regions.  Figure  2(a) 
shows  a  composite  of  two  random  textures  in  three  connected  regions.  In  this  example,  there 
is  no  difference  in  the  mean  gray  level  for  the  two  types  of  textures.  Figure  2(h)  shows  the  Ml. 
segmentation  obtained  by  assigning  pixels  to  black  or  white  levels  (texture'  types  0  and  1,  re¬ 
spectively).  Although  the  outline  of  the  true  regions  is  discernible  to  a  human  observer  the  pro¬ 
cedure  has  in  fact  generated  many  false  regions. 

A  better  procedure  is  to  use  some  form  of  Bayes  estimation.  From  our  arguments  with 

respect  to  the  Ml.  estimate,  we  can  observe  that  a  segmentation  IH  j,  R^,,  ....  R  )  of  the  image 

is  completely  equivalent  to  an  assignment  of  the  pixels  to  0  and  t.  We  shall  refer  to  such  an 

assignment  for  a  point  (ti,  m)  as  its  "state"  s^  ]  )  and  let  that  state  be  stochastically  dependent  on 

the  values  of  states  in  a  surrounding  region  S  .  In  particular,  following  Kaufman  cd  al/  we 

assume  that  the  states  form  a  Markov  chain  with  transition  probabilities  Pi  Is  IS  I.  Bv 

'  1  n,m  n,m' 

combining  Eq.  (1)  and  the  transition  probabilities,  and  applying  Bayes'  rule,  we  can  form  the 
posterior  probability  of  the  state  assignment.  Maximizing  this  expression  leads  to  a  MAI’  esti¬ 
mate  for  the  regions.  By  taking  logarithms  and  eliminating  constant  terms,  we  find  that  an 
equivalent  statement  of  the  MAP  estimation  problem  is  to  minimize  the  expression 
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over  all  choic  es  of  the  s^  .  Specifically,  in  the  two-class  c  ase,  we  must  satisfy  the  conditions 
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The  state  interdependence  imposed  by  the  Markov  model  is  clear  trom  Kq.  (b),  but  it  also  renders 
a  direct  solution  impossible  in  prac  tice.  However,  it  is  possible  to  set  up  an  iterative  proce¬ 
dure  where,  beginning  with  the  Ml,  state  assignments,  one  fixes  the  Sn  and  uses  Icq.  (b)  to  ob¬ 
tain  an  updated  set  of  state  assignments.  The  updated  state's  are  then  used  for  the  S  m  the 

n.m 

next  iteration. 

For  the  results  reported  here,  the  state  was  taken  to  he  dependent  on  value's  in  a  6-  >  *>- pixel 
square  region  centered  around  the  point  (n,  m),  anil  the  transition  probabilities  were  taken  to  he 
proportional  to  the  number  of  l’s  or  0's  present  in  that  region.  Figure  2(c  )  shows  the  results  of 
this  procedure  after  16  iterations.  Since  there  is  an  inherent  ambiguity  in  the  estimation  of  the 
region  boundaries  proportional  to  the  size  of  the  support  of  the  autoregressive  filters,  the  edges 
are  somewhat  rough.  Smoothing  with  a  small-kernel  linear  filter  produces  the  result  shown  in 
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Fig.  2(d).  The  three  regions  have  been  identified  well  except  for  one  small  false  white  area  in 
the  upper  left  corner. 

The  MAP  region  estimation  procedure  was  also  applied  to  the  aerial  photographs  of  terrain 
shown  in  Fig.  3.  The  goal  here  is  to  segment  the  images  into  regions  of  trees  and  fields  without 
further  segmentation  with  respect  to  tree  or  field  characteristics.  Although  both  images  shown 
in  Fig.  4  represent  similar  types  of  data,  there  are  some  significant  differences  between  tin-  im 
ages  particularly  with  respect  to  average  gray  level  of  the  two  classes.  However,  for  purposes 
of  these  experiments,  no  normalizing  procedures  were  applied  to  the  data.  The  filter  coeffi¬ 
cients  for  the  segmentation  and  classification  were  estimated  from  data  in  the  two  white  boxes 
in  Fig.  3(a).  The  MAP  algorithm  was  then  applied  to  the  entire  images.  Figure  4  shows  tin  n 
suits  of  region  estimation  for  the  two  images.  No  smoothing  of  the  edges  was  performed  here. 
Figure  4(a)  shows  generally  solid  identification  of  the  two  region  types.  Figure  4(h)  shows  <  .>n 
sistent  identification  of  the  fields  and  most  of  the  trees  with  some  problems  occurring  m  tree 
areas  that  are  lighter  in  tone  or  show  longer  spatial  autocorrelation  characteristics. 


(a)  (b> 


Fig.  4.  MAP  segmentation  of  Fig.  3  images. 

As  a  follow-on  experiment,  the  upper  left-hand  corner  of  the  image  in  Fig.  tlnl  was  further 
segmented  to  separate  the  larger  trees  from  the  smaller  ones.  In  the  enlarged  .set  inm  :  Imv.i, 
in  Fig.  3(a),  it  is  more  difficult  to  visually  identify  the  image  as  that  of  trees  hut  there  are  now  > 
able  differences  in  spatial  correlations  of  the  textures  within  the  image  that  are  not  unlike  tlm.  i 
seen  in  the  texture  image  of  Fig.  2.  The  segmented  result  shown  in  Fig.  3(b)  tommies  .wth  mo.  i 
observers'  visual  segmentation  of  the  tree  image  and  also  with  the  details  visible  on  the  oririi.nl 
aerial  photograph  prior  to  digitization  (not  shown).  This  suggests  that  terrain  segmentation 
might  well  be  performed  in  a  layered,  hierarchical  manner  starting  with  a  segmentation  into 
gross  categories  and  proceeding  to  increasing  levels  of  detail  as  required. 

2.4  Classification  of  Presegmented  [mages 

In  the  preceding  section,  the  goal  was  to  segment  an  image  into  known  categories  and  the 
segmentation  and  classification  proceeded  concurrently.  In  some  applications,  it  mav  he  require 
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Fig.  5.  Segmentation  of  large  and  small  trees:  (a)  section 
of  tree  and  field  image  and  (b)  MAP  segmentation. 


to  classify  an  image  consisting  of  a  single  type  of  texture  or  terrain  into  one  of  a  set  of  known 
categories.  Such  presegmented  but  unclassified  images  may  arise  as  the  output  of  a  view  of  a 
camera  or  sensor  to  an  area  known  to  be  of  homogeneous  structure.  In  this  situation,  the  clas¬ 
sification  can  be  addressed  in  terms  of  statistical  hypothesis  testing.  Once  again,  the  model 
and  likelihood  functions  described  in  Sec.  2.2  play  a  key  role.  Classification  of  the  image  can  be 
carried  out  as  a  fixed  sample  size  test  involving  the  entire  set  of  pixels  or  as  a  sequential  test 
using  some  contiguous  subset  of  pixels.  The  sequential  test  may  be  useful,  for  example,  if  the 
image  is  being  developed  in  a  scan  mode  and  time  and  resources  for  classification  must  be  care¬ 
fully  allocated.  In  either  case,  the  likelihood  function  is  developed  in  a  recursive  manner.  Addi 
tional  details  of  this  procedure  can  be  found  in  Ref.  2. 

3.  ENHANCEMENT  OF  AERIAL  PHOTOGRAPHY 

3.1  Introduction 

In  this  research  area,  we  have  devoted  our  efforts  to  enhancement  and  restoration  of  aerial 
photographs,  and,  in  particular,  the  sampled  images  recently  received  from  RADC/lRHE.  Our 
studies  can  be  categorized  in  the  following  way: 

(a)  Generalizations  and  applications  of  the  iterative  implementation  of  2-D 
digital  filters.1 

(b)  A  theoretical  study  of  convergence  of  signal  reconstruction  algorithms 
such  as  phase-only  or  magnitude-only  reconstruction,  and  band-limited 
extrapolation. 

(c)  Estimation  of  the  phase  of  the  Fourier  transform  of  an  image  for  image 
restoration. 

In  the  first  area,  we  have  extended  the  original  formulation  of  iterative  filters  so  that  any 
stable  2-D  rational  transfer  function  can  be  implemented  by  an  iterative  filter.  That  is,  we  have 
extended  the  class  of  filters  that  t  an  be  implemented  iteratively.  This  generalization  has  led  to 
some  interesting  filter  structures  for  accelerating  convergence.  In  particular,  trade-offs  be¬ 
tween  filter  complexity  and  speed  of  convergence  have  been  empirically  demonstrated. 

In  a  more  practical  vein,  we  have  applied  iterative  filtering  to  the  problem  of  image  restora¬ 
tion,  and  in  particular,  to  the  problem  of  image  deblurring.  In  addition,  we  also  have  performed 
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some  preliminary  studies  and  experiments  on  the  use  of  iterative  filters  in  cloud  removal,  where 
cloud  cover  is  modeled  by  both  gray-scale  modification  and  blurring  due  to  light  scattering. 

In  the  second  area  of  interest,  we  have  formulated  a  general  convergence  proof,  applicable 

4 

to  a  particular  class  of  signal  reconstruction  algorithms.  Although  our  development  centers  on 
two  specific  examples,  band-limited  extrapolation  and  phase-only  reconstruction,  our  approach 
is  general  and  may  be  applied  to  other  iterative  algorithms  that  satisfy  the  same  assumptions. 

Our  technique  yields  the  first  proof  of  convergence  for  the  phase-only  reconstruction  algorithm 
and  may  be  easily  generalized  to  multi-dimensional  signals.  Furthermore,  the  method  of  proof 
can  be  used  with  nonlinear  constraints  such  as  positivity  and  with  modifications  for  accelerating 
convergence.  Our  studies  also  have  led  to  some  interesting  observations  on  the  relationship  be¬ 
tween  iterative  filters  and  iterative  reconstruction. 

The  final  area  of  investigation  is  phase  estimation,  and  phase  determination  from  intensity 
measurements.  This  work  can  be  classified  as  follows: 

(a)  Determining  the  phase  of  a  complex  wavefront  from  multiple  intensity 
measurements. 

(b)  Estimation  of  the  phase  of  an  image  in  the  presence  of  additive  noise. 

The  first  approach  is  useful  when  a  coherent  optical  imaging  system  yields  a  misfocused 
image.  The  second  area  is  applicable  in  enhancing  images  degraded,  for  example,  by  quantiza¬ 
tion  noise  or  film  graininess. 

3.2  Generalizations  and  Applications  of  Iterative  Filters 

The  original  motivation  for  an  iterative  implementation  of  2-D  digital  filters  was  to  build  up 
a  complicated  impulse  response  from  simple  convolutions  with  impulse  responses  of  finite  ex¬ 
tent.  Such  an  implementation  is  well  suited  to  highly  parallel  array  architectures.  Another  moti¬ 
vation  is  the  implementation  of  non-eausal  rational  transfer  functions.  Image  restoration,  for 
example,  often  requires  a  non-causal  inverse  filter. 

The  iterative  filtering  procedure1  can  be  viewed  as  a  lst-order  feedback  loop,  where  a 
filtered  version  of  the  input  image  is  added  to  a  filtered  version  of  the  result  of  the  previous 
iteration.  A  generalization  of  this  structure  can  be  derived  in  a  way  similar  to  the  original 
derivation. 

3,2.1  A  Generalized  Formulation 

A  2-D  rational  frequency  response  can  be  written  as 

(7) 

where  Afuij.w^)  and  are  trigonometric  polynomials.  In  generalizing  the  original  itera¬ 

tion.  we  define  the  function 

C(iopu>2)  -  1  Mu)p  l®) 

where  is  a  trigonometric  polynomial.  If  we  define  as  the  spectrum  of  the  fil¬ 

ter's  input  signal  xlnj.n^)  and  YicUj.w^)  as  the  spectrum  of  the  filter's  output  signal  yinj.n^), 
then  from  Eqs.  (1)  and  (8)  we  can  derive  an  implicit  relation  for  YftVj.a:^): 

Y(wj,w2)  -  A(o X(cu, ,  +  Clajj.ur,)  Yiouj.a.’^)  .  (9) 
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This  implicit  formula  suggests  the  iterative  formulation: 


Yi(Wj,u)2)  =  X(u)j,w2)  A(a>j,o)2)  X(aij,u2)  +  C(uij,w2)  Y._  j(wj.  u>2)  .  (10) 

The  iterative  solution  will  converge  to  the  desired  filter  output  providing 

|C(o,1,u>2)|  <  1  .  (11) 

In  fact,  unlike  with  the  original  iteration  [where  X(a>1(u>2)  =  1),  it  is  straightforward  to  demon¬ 
strate  that  this  condition  can  always  be  made  to  hold  with  an  appropriate  choice  of  A(a>j,  u>J. 

For  example,  when 

X(aij,a)2)  =aB*(Wj,u>2)  (12) 

where  o  is  a  scalar  such  that 

0  <  a  <  2  max  |  B(olj1,cj2)  |2  (13) 

then,  |C {<*>,,  u>2)|  <  1  for  all  (<i>f,ci>2). 

3.2.2  Methods  of  Accelerating  Convergence 

Equation  (13)  suggests  a  means  of  accelerating  the  convergence  of  the  iteration  given  by 
Eq.  (11).  Since  the  convergence  rate  increases  as  | C'(u)j,  uk)  |  approaches  zero,  we  wish  to  choose 
the  parameters  of  X(oj1,uj2)  so  that 

|C(«1,w2)|  =  |  1  —  Alwj.a^)  B(uj,o>2)  |  »•  0  (14) 


or  equivalently 

A(u>j,  u>2)  »  1/B(o)j,  w2)  .  (15) 

If  we  consider  X(o)j,u>2)  to  be  the  frequency  response  of  a  2-D  filter,  we  can  drive  C(u)j,w2) 
closer  to  zero  by  allowing  the  spatial  extent  of  that  filter's  impulse  response  to  grow.  Thus,  a 
trade-off  exists  between  filter  complexity  and  convergence  rate. 

From  a  slightly  different  point  of  view,  suppose  we  want  to  implement  the  filter  i/B(u>j,u>2) 
and  suppose  A(oij,w2)  is  a  poor  approximation  to  l/B(o>j ,  ca2).  Then,  using  the  iteration  given  by 
Eq.  (11),  we  can  improve  the  result  of  the  filtering  operation  to  achieve  a  better  realization  of 
1/B(c0j ,  co2). 

3.2.3  Application  to  Image  Deblurring 

Let's  now  assume  that  an  image  i(nj,n2)  has  been  blurred  by  a  finite  extent  Gaussian- shaped 
blurring  function  bfn^.n^)  to  give  a  blurred  image  x(nj,n2>.  In  the  frequency  domain,  we  have 


X(Wf,b>2)  =  I(Wj,u)2)  B(W|.w2) 

The  desired  image  is,  of  course, 

Y(wj,o>2)  =  I(ajj,w2)  =  X(a)j,  w2)/B((vj,(a2) 


(17) 


and  thus  the  filter  H(u)j,u)2)  to  be  implemented  is 

H(w1,w2)  =  1/B(o)1,oj2) 


(18) 


An  original  RADC  image  and  its  artificially  blurred  counterpart  are  shown  in  Fig.  6. 
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Fig.  6.  (a)  Original  and  (bj  blurred  RADC  image. 


Since  B(wj,u>2)  takes  on  very  small  values  in  high-frequency  regions,  will  take  on 

exceedingly  large  values  in  these  regions,  and  consequently  will  be  very  sensitive  to  any  noise. 
The  result  of  applying  a  direct  inverse  filter  to  is  shown  in  Fig.  7  along  with  the  blurred 

image.  Note  the  severe  degradation  due  to  high-frequency  sensitivity. 


TTmoii- »] 
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Fig.  7.  (a)  Blurred  RADC  image  and  (b)  reconstruction  by  direct 

inverse  filtering. 


The  iterative  approach,  on  the  other  hand,  slowly  builds  up  the  desired  inverse  filter  and 
seems  to  avoid  the  sensitivity  problem  inherent  in  the  direct  approach.  Figure  8(a)  shows  a 
blurred  image  and  (b)  illustrates  the  restored  image  after  20  iterations  with  set  equal 

to  a  constant.  When  X(u.'(,u>2)  is  allowed  to  vary  (a  filter  of  length  equal  to  that  of  B(iuj,u'2)  was 
designed),  fewer  iterations  were  required  to  obtain  a  slightly  better  restoration.  Figure  9(a) 
depicts  the  blurred  imane  and  (b)  the  restored  image  after  only  10  iterations  when  the  variable 
A(u.j,u)2)  was  applied.  This  example  demonstrates  the  trade-off  between  the  spatial  extent  (or 
complexity)  of  and  the  convergence  rate. 
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Fig.  8.  (a)  Blurred  HADC  image  and  (b)  reeonstruetion  by  iterative 

filtering  with  fixed  A  and  20  iterations. 


(o)  (b) 


Fig.  9.  (a)  Blurred  HADC  image  and  (b)  reeonstruetion  by  iterative 

filtering  with  variable  A  and  10  iterations. 


! 

3.2.4  Extensions 

The  theoretical  and  empirical  results  of  the  previous  subsections  have  led  to  a  number  of 
interesting  possible  extensions.  First,  under  certain  conditions,  the  generalized  formulation 
can  be  shown  to  yield  the  phase  of  the  desired  (filtered)  image  after  only  one  iteration.  Since, 
for  a  certain  class  of  images,  all  the  information  lies  in  the  phase,  we  might  conclude  that  after 
only  one  iteration,  we  have  restored  the  desired  information.  This  phase  information,  therefore, 
may  be  useful  in  increasing  the  convergence  rate  or  perhaps  in  reducing  the  noise  buildup  as  the 
iteration  proceeds. 

Another  important  insight  we  have  recently  made  is  that  the  iteration  of  Eq.  (11)  can  be 
slightly  modified  to  lake  on  the  appearance  of  the  more  traditional  iterative  signal  reconstruction 
•  algorithms  such  as  the  Gerchberg5  and  Van  Cittert^1  algorithms.  As  a  result,  it  is  possible  to 

view  an  iterative  filter  as  an  iterative  reconstruction  algorithm  so  that  procedures  used  within 
these  more  traditional  approaches  can  be  applied  to  our  problem.  For  example,  we  are  consid- 

'  ering  imposing  finite  extent  and  positivity  constraints  on  the  output  of  each  iteration. 

i  i 
S  ; 

[  |  3.3  Convergence  of  Iterative  Signal  Reconstruction  Algorithms 

| ! 

'  Imposing  finite  extent  and  positivity  constraints,  as  described  in  the  previous  subsection, 

leads  into  issues  of  the  convergence  of  a  modified  version  of  the  iteration  given  by  Eq.  (11). 

More  generally,  there  exist  many  iterative  procedures  for  signal  reconstruction  which  impose 
such  constraints  in  the  space  domain  and  also  makes  use  of  known  information  in  the  frequency 
domain.  For  example,  often  only  the  low-pass  region  of  the  spectrum  of  an  image  is  given.  In 
this  case,  we  would  attempt  to  retrieve  resolution  by  iteratively  incorporating  the  low-pass  in¬ 
formation  in  the  frequency  domain,  and  the  known  extent  and  positivity  in  the  space  domain. 
Furthermore,  if  in  addition  to  being  low-pass  filtered,  the  data  are  also  degraded  by  a  known 

zero-phase  blurring  function,  only  the  phase  of  the  low-pass  region  may  be  accurately  preserved. 
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An  iteration  had  been  developed  for  recovering  the  desired  image  from  this  phase  information.  ’ 

■I 

More  recently,  we  have  developed  a  proof  of  convergence  for  this  phase-only  iteration.  The 
generality  of  the  proof  allows  for  nonlinear  constraints  such  as  positivity  and  minimum  or  max¬ 
imum  value  constraints,  and  also  encompasses  a  large  class  of  iterative  procedures  which  are 
based  on  a  nonexpansive  mapping. 

3.3.1  Nonexpansive  Reconstruction  Algorithms 

The  convergence  proof  relies  on  the  property  that  an  iterative  reconstruction  procedure  is 
based  on  a  nonexpansive  mapping.  Heuristically,  a  nonexpansive  mapping  defined  on  a  set  of 
images  has  the  characteristic  that  when  any  two  images,  x  and  y,  are  operated  on  by  the  map¬ 
ping,  the  resulting  images  are  "closer"  to  one  another  in  a  mean-squared  sense,  l.etting  d(x,y) 
denote  the  mean-squared  distance  between  two  images,  we  say  a  mapping  F  is  nonexpansive  if 

d(Fx,  Fy)  <  d(x,y)  .  (19) 

(Fx  denotes  the  mapping  F  applied  to  the  image  x.) 

The  class  of  algorifhms  we  have  dealt  with  are  strictly  nonexpansive;  that  is, 

d(Fx,  Fy)  <  d(x,  y)  when  x  ^  y 

This  strict  nonexpansiveness  property  has  been  proven  for  both  the  phase-only  reconstruction 

4 

iteration  and  also  the  band-limited  extrapolation  iteration.  More  generally,  it  is  of  major 
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importance  in  proving  convergence  of  a  class  of  iterative  algorithms  of  the  form 


xk  +  1  =  Fxk  .  U«» 

Moreover,  it  marks  a  significant  difference  from  the  contraction  property  which  has  been  widely 
used  in  the  literature  for  proving  the  convergence  of  iterative  procedures  based  on  contraction 
mappings. 

3,3.2  Convergence  of  the  Under- Helaxed  Iteration 

One  important  consequence  of  our  approach  is  that  convergence  can  be  demonstrated  for  a 
modified  form  o  Eq.  (20)  (Refs.  4  and  9).  For  example,  the  phase-only  iteration  can  be  modified 
to  accelerate  convergence  by  forming  the  "under- relaxed"  version  of  the  iteration  given  by 

xk+l  =  11 -X>  *k  +XF\  ■  U1> 

with  appropriate  choice  of  X.  Eq.  (21)  converges  faster  than  the  iteration  (Eq.  (20) |.  The  gener¬ 
alized  iteration  with  variable  X  also  can  be  put  in  the  form  of  Eq.  (21).  As  a  result,  the  conver¬ 
gence  of  certain  accelerated  versions  of  the  iterative  filter  implementation  may  be  proved. 

3.4  phase  Estimation  for  Image  Restoration 

The  final  topic  to  be  discussed  involves  techniques  of  determining  phase  information  from 
intensity  measurements.  We  have  derived  a  recursive  procedure  for  generating  phase  informa¬ 
tion  from  intensity  measurements  in  more  than  one  plane.  The  derivation  is  given  in  1-D,  but 
it  is  easily  extended  to  2-D.  This  solution  has  potential  application  to  focusing  a  misfocused 
image  formed  by  a  coherent  imaging  system. 

Suppose  a  complex  signal  x(n)  is  the  input  of  a  linear  space-invariant  filter  with  impulse 
response  h(n).  Then 

y(n)  =  ^  h(n  -  k)  x(k)  (22) 

k 

where  y(n)  is  the  filter  output.  When  the  magnitude  functions  |y(n)|  and  |x(n)|  arc  known,  a  re¬ 
cursive  algorithm  can  be  derived  for  generating  all  possible  sequences  x(n)  and  thus  all  possible 
phase  functions  associated  with  lx(n)  |. 

The  recursive  solution  is  described  in  the  following  way.  Letting  r  and  i  denote  "real  and 
imaginary  part  of,"  respectively,  we  must  solve  the  set  of  equations  for  x  (n)  and  x.(n)  given  by 


x^ln)  dr(n)  -  xt(n)  d.(n)  =  e(n) 

(23  a) 

|x(n)  \l  -  x^(nj  +  x^(n' 

(2  3b) 

where  x(n)  =  x^n)  +  jXj(n),  and  where  d  (n),  d.(n),  and  c(n|  are  functions  of  h(n),  Ixtn)',,  ly(n)|, 
and  the  previous  values  of  x(n),  i.e.,  x(m),  m  -  0,  t...n  -  t.  It  is  straightforward  to  see  from 
Eqs.  (23a)  and  (23b)  that  there  exist  two  possible  solutions  for  x(n).  We  choose  one  of  the  two 
possible  values  of  x(n)  and  proceed  to  compute  x(n  +  t).  Different  selections  of  x(n)  at  each  step 
will  result  in  different  estimates  of  x(n),  but  sequences  which  are  compatible  with  the  known 
magnitude  functions  !x(n)|  and  |y(n)  |. 


One  approach  to  removing  this  two-solution  ambiguity  is  to  assume  the  presence  of  a  second 
filter  acting  upon  x(n),  and  that  the  magnitude  of  its  output  is  known.  With  this  second  filter,  a 
necessary  and  sufficient  condition  on  the  two  filters  was  derived  for  generating  a  unique  solution 
for  x(n). 

One  potential  application  of  the  recursive  solution  is  the  reconstruction  of  a  complex  propa¬ 
gating  coherent  wavefront  from  intensity  measurements  in  two  or  more  planes.  The  2-0  filter 
in  this  case  represents  the  propagation  phenomenon  between  planes. 

The  recursion  can  be  used  to  generate  the  complex  waveform  along  one  plane.  The  2-D  fil¬ 
ter  can  then  be  used  to  "propagate"  coherent  light  to  any  other  plane  in  space.  For  example,  if 
an  image  along  a  particular  plane  is  out  of  focus,  the  filter  simulation  can  theoretically  bring 
the  image  into  focus  at  some  other  plane. 
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